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Abstract. 

We present a study of binary mixtures of Bose-Einstein condensates confined in a 
Q ' double- well potential within the framework of the mean field Gross-Pitaevskii equation. 

O/ We reexamine both the single component and the binary mixture cases for such a 

potential, and we investigate in which situations a simpler two-mode approach leads 
to an accurate description of their dynamics. We also estimate the validity of the 
most usual dimensionality reductions used to solve the Gross-Pitaevskii equations. 
r — ^ I To this end, we compare both the semi-analytical two-mode approaches and the 

O^ ' numerical simulations of the ID reductions with the full 3D numerical solutions of 

the Gross-Pitaevskii equation. Our analysis provides a guide to clarify the validity 
^^ ' of several simplified models that describe mean field non-linear dynamics, using an 

o 



experimentally feasible binary mixture of an i^ = 1 spinor condensate with two of its 
Zeeman manifolds populated, m — ±1. 



PACS numbers: 03.75.Mn 03.75.Lm 03.75.Kk 74.50.+r 
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1. Introduction 

The phase coherence of a Bose-Einstein Condensate (BEC) is an important and 
characteristic property of ultracold bosonic gases that leads to fascinating macroscopic 
phenomena such as interference effects or Josephson-type oscillations. Two condensates 
trapped in a double-well potential exhibit interference fringes when the barrier is released 
and the two expanding condensates, with a well-defined quantum phase, overlap. 
Instead, if the barrier is not switched-off and is large enough to ensure a weak link 
between both condensates in each side of the trap, the quantum phase difference will 
drive Josephson-like effects, which consist on fast oscillating tunneling, much faster than 
the single particle tunneling, of atoms through the potential barrier [H [2]. 

The first evidence of the phase coherence of a BEC was obtained in early interference 
experiments [3j where clean interference patterns appeared in the overlapping region of 
two expanding condensates. It has been only recently that a clear evidence of a bosonic 
Josephson junction in a weakly linked scalar BEC Ij has been experimentally reported 
by the group of M. Oberthaler in Heidelberg |3]. In this experiment, two condensates 
are confined in a double-well potential with an initial population imbalance between 
both sides which triggers the Josephson oscillations. The tunneling of particles leads 
to a coupled dynamical evolution of the two conjugate variables, the phase difference 
between the two weakly linked condensates and their population imbalance. In spite 
of the system being very dilute, the inter-species interaction plays a crucial role in the 
Josephson dynamics, leading to new regimes beyond the standard Josephson effect, e.g. 
macroscopic quantum self trapping (MQST). 

The Gross-Pitaevskii (GP) mean field theory provides a natural framework 
for investigating Josephson dynamics in weakly interacting systems at very low 
temperature. Josephson oscillations in scalar Bose-Einstein condensates have been 
theoretically studied by using different techniques |5l|6|[71[8|[9|[Tnt[TTt[T2l[T8|[Til 
[T5l [T6| Wl\ [T8| [T9] . As expected, the full three dimensional time-dependent Gross- 
Pitaevskii equation (GP3D ) provides an excellent agreement with the experimental 
data [9l[T0l[20]. However, since 3D dynamics need in general rather involved calculations, 
one can benefit from the fact that the barrier is created along one direction and the 
tunneling of particles is mainly one dimensional (ID) to investigate the Josephson 
dynamics by means of effective ID GP-like equations. Among these reduced GP 
equations, the non-polynomial nonlinear Schrodinger equation (NPSE ) proposed in 
Ref. [21] has provided the best agreement with the experimental results in scalar 
condensates [20], whereas another effective ID Gross-Pitaevskii equation ( GPID ) 
fail to describe the dynamics for large number of trapped atoms in the same trapping 
conditions as in the Heidelberg's experiment [TOl |20] . 

Interactions are important to understand the different regimes of the tunneling 
dynamics. Therefore, multi-component BECs in double-well potentials offer an 
interesting extension to study phenomena related to phase coherence. In particular 

X The notation "scalar BEC" is used as equivalent to "single component BEC" in this article. 
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the Josephson dynamics will become richer due to the interplay between intra- and 
inter-species interactions. 

Josephson oscillations in binary mixtures confined in double-well potentials have 
been addressed in a number of recent articles. The case of two-component BECs 
with density- density interactions has been studied within two-mode approaches in 
Refs. [22l[23l[2a[25l|26l|271|28l|29]. Refs. [261 |29] go one step further and also 
consider GPID simulations. Spin-dependent interactions have been addressed in 
Refs. [301 SI]- Josephson dynamics in spinor condensates confined in double- wells, 
characterized by an exchange of population between different Zeeman components, 
has also been investigated in Refs. [321 |33]. In Refs. [26l |30] the interest of studying 
Josephson dynamics in binary mixtures has been emphasized as it can give access to 
information of the different scattering lengths present in the system. 

Recently, the equations of the tunneling dynamics in a binary mixture within the 
two- mode approximation to the GP equations have been derived in Ref. [25]. However, 
the authors have not compared their two-mode analysis to direct numerical resolutions of 
the GP equation, and have also not provided microscopic values to the parameters of the 
two-mode equations. Their main result is the description of a symmetry breaking pattern 
occurring when the inter- and intra-species interactions differ substantially. In Ref. [26] 
a comparison of the standard two-mode approach and the coupled GPID equations for 
the mixture has been presented for one specific double-well potential which allows an 
analytical treatment. 

For single component BECs, it has been already studied the range of validity 
of the different approximations to the Josephson dynamics by comparing with 
GP3D calculations and with the experimental results. However, no comparison with 
the full GP3D dynamics has been yet performed for a binary mixture in a double-well 
potential. 

The aim of this paper is to investigate systematically the tunneling dynamics of 
a binary mixture of BECs trapped in a double-well potential, as well as the validity 
of the different mean-field approximations. We consider a mixture of two components 
obtained by populating two Zeeman states of an F = 1 *^Rb condensate confined in 
the same double- well potential as in the experiments [4]. This system corresponds to a 
natural extension of the experimental work of Ref. [1], where only one of the Zeeman 
components was populated. 

We provide a general overview of the different techniques used to investigate 
Josephson dynamics within the two-mode model (standard and improved two-mode) 
and within the Gross-Pitaevskii framework (one-dimensional reductions of the GP 
equation, GPID and NPSE). To this end, we solve the full 3D time-dependent GP 
equation for the mixture as a reference to assess and analyze the validity of the previous 
approximations. 

The paper is organized as follows. The general framework of the coupled Gross- 
Pitaevskii equations for a binary mixture is presented in Sec. O In Sec. [3l we derive 
analytic two-mode models both for single and two component systems. First we 
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recall the standard two-mode model (S2M). Then we derive the equations of the 
improved two-mode model (I2M ) for a binary mixture, generalizing the work for a 
single component BEC performed in Ref . ^ . We discuss the stability of the dynamical 
equations and look for the stationary points for a binary mixture. In Sec. HJ we 
analyze the different one-dimensional reductions of the GP3D equations for the mixture: 
GPID and NPSE . In Sec. O we revisit the dynamics of a single component condensate 
in a double- well potential with the same parameters as in the experiment |1]. The 
tunnehng dynamics in two-component systems is accurately discussed in Sec. El We 
obtain the dynamics by solving the coupled GP3D equations for the mixture and show 
that for certain conditions there exists a good agreement between I2M and GP3D , as 
well as for NPSE and GP3D . The range of validity of the two-mode models is explored, 
paying special attention to situations which fall beyond the two-mode approximation. 
Finally, we discuss cases that present characteristic features arising from the mixture, 
with no analog in the tunneling dynamics of a single component BEC. Conclusions are 
given in Sec. [71 

2. Mean field approach: Gross-Pitaevskii equations 

We consider a binary mixture of weakly interacting atoms at zero temperature, confined 
by the same double- well potential, V^(r). For dilute systems with sufficiently large 
number of particles, the Gross-Pitaevskii equation provides a suitable framework to 
study the dynamics. In the mean field approximation, each condensate is described 
by the corresponding wave function \E'j(r;t), with i = a,b denoting each of the two 
components of the binary mixture. To avoid any misunderstanding let us remind the 
reader that we are describing two different kind of atoms, a and b, which evolve on a 
double- well external potential. In most situations, the system will behave as if there were 
four weakly linked Bose-Einstein condensates, two per each component of the binary 
mixture per each side of the potential barrier. The mean field description will refiect 
this feature by the homogeneous quantum phase of \E'j(r; t) at each side of the potential 
barrier, as will be discussed in great detail in the following sections. 

The dynamical evolution of the two wave functions can be obtained by solving the 
two coupled GP equations: 



ih 
ih 



d^a{r;t) 

dt 
d ^b{r;t) 

dt 



h^ 



2ma 
2mh 



V' + V{r) + g,aNa\^,{r;t)\^ + g,,N,\^il,{r;t)\^ 



*a(r;t) 



^fe(r;t) 



For each component, the condensate wave function \I'i(r;t) is normalized to 1, rrii is 
the atomic mass, and ga = Anh'^ai/rrii is the effective atomic interaction between atoms 
of the same species, with Oj the corresponding s-wave scattering length. The coupling 
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between both components is governed by the inter-species interaction Qah = Qba-, which 
depends on the specific nature of the binary mixture. The total number of atoms in the 
mixture is A^ = A'^^ + A*";,. 

There are many experimental possibilities to study the dynamics of binary mixtures 
of BECs. We will restrict our study to one of them, which is experimentally feasible. 
We will consider binary mixtures made of F = 1 ^^Rb atoms populating the m = ±1 
Zeeman sublevels [31] . This implementation greatly simplifies the dynamics as the inter- 
and intra-species couplings are very similar in magnitude. Of course this choice limits 
the phenomena which can be observed, e.g. the interesting symmetry breaking pattern 
discussed in Ref. [25], which relies on the inter-species coupling being larger than the 
intra-species one, will not take place, see Sec. 16.51 

On the other hand its simplicity allows to discuss in detail the different approaches 
taken in the literature, e.g. two-mode models of the GP equations, one dimensional 
reductions, etc. As occurred in the scalar case, the dynamical features contained in 
Eqs. ([T]) can to a large extent be described by a simplified two-mode model for each 
component. In the next section we follow Refs. [5], [9] and [25] and derive two- mode 
expressions for the scalar and binary case. The usual assumption of neglecting the 
overlaps involving the right and left modes gives rise to the so called standard two-mode 
(S2M ) equations, while retaining them one also gets a closed system of equations, the 
improved two- mode (I2M). Both the S2Mand I2Mare also derived for the binary 
mixture case. 

3. Two-mode approaches 

The two-mode approximation allows to study the dynamics of weakly linked Bose- 
Einstein condensates, without solving the full GP3D nor reducing the dimensionality 
of the GP equation [H [9]. Depending on the specific double- well potential, e.g. on 
the energy gap between the first two levels of the single particle Hamiltonian and the 
next two, it can provide an excellent description of the full GP solution. The relevant 
physical quantity is the ratio between the energy gap between the ground state and first 
excited state of the double-well potential, 5Eq.i = Ei — Eq^ and the energy difference 
of the ground state and the second excited state, 5£'o;2 = E2 — Eq. The smaller the 
ratio 6Eq-i/6Eq-2 the more accurate the two-mode approach. The two mode description 
characterizes the dynamics of the scalar condensate in a double-well potential with only 
two variables: the relative population and the phase difference between the left and 
right side of the potential barrier. 

§ Note that this is zero if the barrier is infinitely high. 
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Figure 1. 3D depictions of the density, p{x,y) = J dz\^{x,y, z)\'^ in {fim)~'^, of the 
(a) ground state, (6) first excited state, (c) left mode and (d) right mode obtained by 
performing an imaginary time GP3D calculation with the same conditions as in the 
experimental setup of Ref . [1] . 



3.1. Standard two-mode model for the single component case 

The GP equation for the scalar case corresponds to a particular limit of the GP equations 
for the binary mixture, Eqs. ([T]), 



ih 



~~di 



2m 



^(r;t). 



(2) 



We will make use of the following notation, Hq = — (/i^/2m)V^ + V{r), and -f^gf^I/] = 
gN\'^{r;t)\'^. Let us recall the two-mode approximation for a single component 
condensate in a double-well potential. We consider A^ interacting atoms with atomic 
mass m, and coupling constant g, trapped in a symmetric double-well potential V{r). 
When both sides of the potential barrier are weakly linked, the total wave function 
can be approximately written as a superposition of two time-independent spatial wave 
functions $2,(ij)(r) mostly localized at the left (right) side of the trap: 

^(r; t) = ^L(t)$L(r) + ^ij(t)$i,(r) . (3) 

The left and right modes, can be expressed as linear combinations of the ground (+) 
and the first excited (— ) states of the double- well potential including the interaction 
term. They satisfy, {Hq + Hg[^±])^± = /i±$±, and the left/right modes can be written 
as [9]: 

$+(r) + $_(r) , _ $4-fr)-$_fr' 



*L(r) 



V2 



^R{r) 



V2 



(4) 
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We observe that in a symmetric double- well, $-|- have a well defined parity: $±(r) = 
±$-|-(— r) , and therefore ($j$j) = 5.y with i,j = +,—. Since they are stationary 
solutions of the GP equation, ^± are real functions, and so are the left and right 
modes ^l{r)- The integrated density in the z-direction, p{x,y) = J dz\'i/{x,y,z)\'^, 
associated to the ground ($+), and first excited ($-) states are depicted in Fig. [1] 
together with the densities associated to the left and right modes. The plots correspond 
to the experimental set up of Ref . |1] . 

From the phase coherence properties of a BEG, one can assume that the wave 
function in each side of the trap has a well defined quantum phase 0j(t), which is 
independent of the position but changes during the time evolution. We can write. 



M/,(t) = ^iV,(t)e^<^^W, (5) 

where iVi(_R)(t) corresponds to the number of atoms on the left (right) side of the trap, 
and the total number of atoms is A^ = NL{t) + Nji{t). The weak link condition is fulfilled 
if (/i_ - fi+) « (l/2)(/i+ + /i_). 

As a first step, we consider the so-called standard two-mode approximation (S2M ), 
which neglects a certain set of overlapping integrals involving mixed products of $l and 
$/j. This approximation yields essentially the correct qualitative results in the scalar 
condensate although it may lead to incorrect quantitative predictions depending on the 
specific barrier properties [H [9] . 

Inserting the two-mode ansatz ([3]) in the GP equation for a single component 
condensate ([2]) and neglecting terms involving mixed products of $l and $/? of 
order larger than one, yields into a system of equations for the two localized modes 
which can be written in terms of two dynamical variables: the population imbalance 
z(t) = [NL{t) — Nji{t)]/N and the phase difference S(j){t) = (puit) — (piif) between each 
side of the barrier: 



z{t) = -ujR^l-z'^{t) sin5(t){t) (6) 

5<p{t) = unAE + u /^^f'^ Nzit) + uJn^=^B== cos 6<P{t) , 

where ur = 2K/h is the Rabi frequency and 



2K AK 

-.2 



E'uR) = I dv 



I dr[^\V^Hn){r)f + ^l^n){r)V{r) 



K = - dr 



2 



|- V^lW ■ V^iilr) + <l>i(r) V{r) $^(r) 
2m 

UL(R)=g jdv^\^^^{v). (7) 

For a symmetric double- well, E^ = E^ and Ul = Ur = U ■, therefore AE = 0. Moreover, 
the Rabi frequency only appears as a scale in the problem and thus can be absorbed in 
the time by rescaling t — )■ u^t. Then, together with the definition A = NU/{hwR), we 
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z{t) = - VT^^2(^ sin(50(t) (8) 

zit) 
5<P{t) = Kz{t) + ^' cos(50(t) . 

^l-z^{t) 

Note that, A > and A < correspond to repulsive and attractive atom-atom 
interactions, respectively. There are different regimes depending on the initial values of 
the population imbalance and phase difference, z{Q) and (5(/)(0), Sec. 13.31 
From the energy functional of the GP equation (|2]): 

-2m 

and using the two- mode ansatz ([3l), we can define the conserved energy per particle of 
the system as: 

(10) 

where C is a rescaling constant. If we consider again a symmetric double-well we have: 
A 



E[^{v-t)]= /rfr|£-|Vv^(r;t)| + F(r)|vl/(r; t)|V ||vl/(r; t)r | (9) 



H = -z\t) - ^/ 1 - z^{t) cos 5^{t) . (11) 

Note that the equations of motion (JHl) can be written in the Hamiltonian form: 

being z and 6(f) canonical conjugate variables. 

3.2. Improved two-mode model for the single component case 

Ananikian and Bergeman [H] noticed that for a symmetric double-well there was no need 
to neglect any of the overlapping integrals to obtain a closed set of equations relating z 
and 6(f). 

Thus, remaining in the two-mode approximation but retaining all the overlaps it is 
straightforward to write down the following set of equations (cf. Eqs. (22) in Ref. [9]), 
called the "improved two-mode" (I2M ) equations, 

z{t) = - By/1 - z^{t) sm6(f){t) + C{1 - z'^it)) sm26(f){t) (13) 

6(f^{t) = Az{t) + ^^^^ cos(50(t) - Cz{t) cos2(50(t) . 

y/l-Z^{t) 

Defining, jij = g Jdr $2(r)$2(r) ,i,j = +, -, we have, A = N (107+_ - 7++ - 7__) /4, 
B = 2K + Ar(7__ - 7++)/2, and C = gN Jdr $|(r)$|(r). 

As discussed in detail in Ref. [9] , the physics arising from the I2M is similar to 
the one present in the S2M . The I2M , however, is in much better agreement with the 
GP3D for a broader set of double-well potentials, as we will see in Sec. [5l In particular, 
for the double-well considered in the experimental setup of the Heidelberg group [1] 
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Figure 2. Energy surface, Eq. (fTT|) . for A = 2.5. The lines on tlie surface correspond 
to possible trajectories of the system. 



the S2M (with the corresponding microscopic parameters, K and t/, computed from 
the GP3D with the experimental 3D potential and experimental coupling, g) does not 
correctly predict the physics of the experiment, mainly due to the importance of the 
neglected overlaps, and not to a dynamics far from a two-mode one. 



3.3. Regimes for the single component case 

3.3.1. Stability analysis 

In this section we use the S2M to analyze the stability of the single component 
system, and focus on the case of repulsive interactions A > 0. Using the Hamiltonian 



( TTTl) and the equations of motion ( 112 
solving the equations: 



the stationary points 



dH 
dz 







dH 



85(1) 



0. 



^) can be found by 



(14) 



To asses the stability of these points, we need to study the Hessian matrix of the system, 
which for the possible values of the phase difference, 50° = or vr, is always diagonal and 
its eigenvalues are d'^H\zO^S(t>o and dj^H\zO^S(j}0- Depending on the sign of these eigenvalues 
the stationary points will be maxima, saddle points or minima. The stationary points 
and their stability are summarized in Table [1] 

The evolution of the system can be represented on a z — 6(f) plane, where the system 
follows trajectories with constant energy, H, see curves in Fig. [2l Note that oscillations 
around a stationary point, closed curves, occur only if the central point is either a 
maximum or a minimum of the energy, but not a saddle point. As we will see in the 
following sections, these orbits will give rise to the Josephson oscillations and to the 
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Figure 3. Different regimes for a set of initial conditions, imbalance z(0) in the y-axis 
and phase difference <50(O) in the x-axis. The upper panels correspond to repulsive 
interactions while the lower ones to attractive interactions. The values of |A| are 0.5, 
1.5, and 5. for the left, middle and right panels respectively. Grey regions correspond to 
Josephson oscillations, blue regions to 7r-modes (upper panels) and zero-modes (lower 
panels), and red regions to running phase modes. 



zero- and 7r-modes. 

3.3.2. Symmetry between attractive and repulsive interactions 

The stability analysis has been presented only for repulsive interactions, but from 
the system (El) we can see that if we change the interactions, A — y —A, we recover the 
same system of equations if 6(f) — > vr — 6(f): 






- Vl-z2(t) sin (vr - (50(t)) 
z{t) 



(15) 



cos(7r — 6(f){t)) 



which means that the dynamics of the system and the different regimes are the same 
for both interactions, with a phase-shift of vr. This can be seen in Fig. [3l that shows 





(^0, <50O) 


stationary 


minimum 


saddle 


maximum 




(0,0) 


VA 


VA 


— 


— 




(0,7r) 


VA 


— 


A>1 


A< 1 


(±v/i- 


- l/A^ vr) 


A> 1 


— 


— 


A> 1 



Table 1. Stationary points of the system for repulsive interactions, A > 0, and their 
stability. 
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Figure 4. z — </> representation of different constant energy trajectories for three values 
of A: 0.5 (a), 1.5 (b), and 5 (c). Solid-black lines correspond to Josephson oscillations, 
dotted-blue to Tr-modes, and dashed-red lines to running phase modes. 



the behavior of the system for a given set of initial conditions. The upper panels are 
for repulsive interactions A > and the lower ones for attractive interactions A < 0. 
The grey regions correspond to Josephson oscillations, the blue regions to zero- and 
TT-niodes, and the red regions to running phase modes, as it will be seen in the following 
sections. 



3.3.3. Josephson dynamics 

This regime is characterized by fast oscillating tunneling of population across 
the potential barrier. Plotted in a z — 50 map, the system evolves following closed 
trajectories around a minimum or a maximum {z^ = 0, 50°) configuration, with a zero 
time-average of the population imbalance, < z >t= 0. The stability analysis shows 
that for A > —1, which corresponds to repulsive or slightly attractive interactions, 
the stationary point {z^ = 0, (50° = 0) is a minimum permitting Josephson oscillations 
around it. Analogously, for A < 1, either attractive or slightly repulsive interactions, 
the stationary point (z° = 0, 50° = it) becomes a maximum, and therefore also allows 
for closed orbits around it. For |A| > 1, there are Josephson oscillations around only 
one point: {z^ = 0, 50° = 0), or (z° = 0, 50° = vr). However, in the region of weak 
interaction, |A| < 1, the oscillations around both points are allowed. 

In panel (a) of Fig. HJ A = 0.5, the black closed orbits around 50° = or around 
50° = 71 correspond to Josephson dynamics around these points. In panel (b) however, 
as A = 1.5 > 1, only the origin can give rise to Josephson oscillations, so the closed 
orbits around (2° = 0, 50° = vr) disappear. 

It is also interesting to study the behavior of the system for small oscillations around 
these two stationary points of zero imbalance, smallest orbits in Fig. H] (a). In this 
limit, the system ([8]) can be linearized giving the dynamical equation: z{t) = —z{t){l + 
A cos 50°) with cos 50° = ±1. The population imbalance performs sinusoidal oscillations 
with a frequency uj = ur^^I + A cos 50°, independent of the initial conditions. Note 
that this frequency only exists when these points are either maxima or minima. The 
phase difference oscillates with the same frequency but with a phase-shift of 7r/2 with 
respect to the imbalance. If the initial population imbalance increases, the dynamics 
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of the system changes substantially to non- sinusoidal oscillations, and the frequency 
becomes dependent on the initial conditions. 

3.3.4- Macroscopic quantum self trapping 

In the case of repulsive interactions, we have seen that for A > 1, the stationary 
point {z^ = 0, 6(f)^ = 7t) becomes a saddle point and there appear two maxima, 
[z^ = ±a/1 — l/A^, 50° = tt). A similar behavior is found for attractive interactions. 
These stationary points allow for oscillations around them with < z >t^ 0. In fact, in 
this regime, the imbalance has the same sign during the evolution, and therefore one of 
the wells is always overpopulated. 

This regime is called macroscopic quantum self trapping, as the tunneling is strongly 
suppressed and the particles remain mostly trapped in one of the wells. This is a 
phenomenon arising from the atom-atom interaction, which appears as a non-linearity 
in the Gross-Pitaevskii equation. 

The critical condition for the existence of the MQST regime can be found by 
imposing that the system remains on one side of the trap [5]. For a given set of 
initial conditions, z{0) ^ and 50(0), the system will remain trapped if. 



A > 2 ^ ' .^ for A > 1 



Vl- 


-z(O)2cos[50(O)] + l 


Z(0)2 


Vl- 


-^(O)2cos[(50(O)]-l 



where the limits of the interaction parameter are due to the fact that only when |A| > 1 
the (z° 7^ 0, (50°) stationary points exist. 

In this regime however, there are two different kind of MQST depending on whether 
the phase difference evolves bounded, giving the so-called zero- and vr-modes, or whether 
it evolves unbounded, increasing (or decreasing) always in time, giving rise to the 
running phase modes. 

For values of the interaction parameter of 1 < |A| < 2 the only MQST regime that 
one can have is the zero-mode for attractive interactions and the vr-mode for repulsive 
interactions (which are plotted in blue dotted lines in panel (b) of Fig. H]). In these 
regimes the phase difference evolves bounded around 50 = and 50 = vr, respectively. 

On the other hand, for values of |A| > 2 one can have both classes of MQST. In 
general however, for a given set of initial conditions, the system will evolve following 
a running phase mode (dashed- red lines of panel (c) of Fig. H]), because the values of 



z^ = ±a/1 — l/A^, that allow closed orbits, are very close to 1 (see the small vr-modes 
of panel (c) in blue dotted lines). 

In panel (c), one can see that the broadest closed orbit around {z^ ^ 0, 50° = vr), 
for A > 2, is the one that goes through {z = ±1, 50 = 0). Notice that an orbit that 
crosses the 50 = axis in any other point, z ^ ±1, would correspond to a running phase 
mode. The case of attractive interactions can be understood by taking into account the 
phase-shift of tt in 50. The latter can be used to find the condition to have bounded 
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or running phase difference modes. For a given set of initial conditions (-2(0), (50(0)) 
fulfilling the self-trapping condition flT6|) . the system will evolve in a bounded phase 
mode only if: 

|A| < ^p% . (17) 

Moreover, in a zero- or a 7r-mode MQST, we can study small oscillations around 
the corresponding minima or maxima, z{t) = z^ + 5z and 5(j){t) = 5(f)^ + S(j){t), so the 
linearized system ([8]) becomes: 



Sz{t) = -6z{t) 



1 + A cos ( 



which gives a sinusoidal behavior with a frequency: 



:i8) 



1 - 2(z'^V 
u = URjl + Acos6r _ J . (19) 

3.4- Standard two-mode model for the binary mixture 

Let us recall the two-mode approximation for weakly linked binary mixtures [22| [23l [2l| 
[25] . The total wave function of each component is written as a superposition of two 
time-independent spatial wave functions localized in each well: 

vl/,(r; t) = vI/,x(t)$,L(r) + ^,RmM^) , (20) 

with ($jQ,|$j^) = 6ij6ai3, i,j = a,b and a,f3 = L,R. For a given component, the 
condensates in each side of the trap are weakly linked. Then, as in the scalar case, one 
can assume that the wave function in each side of the trap has a well defined quantum 
phase (j)j^a(t), which is independent of the position but which changes during the time 
evolution. Thus, 



^,At) = sjNUty'^'^"^'^ . (21) 

Nj,L(R)(^) corresponds to the population of the j-component on the left (right) side of 
the trap, with Nj = Nj^L{t)-\-Nj^ii{t). Inserting the two-mode ansatz fl20|) in the coupled 
GP equations for the mixture ([1]), retaining up to first order crossed terms yields the 
following system of coupled equations: 

2K„ 



Zait) = -—^^l-zlit)sm6Ut) 

h 

, ^Kg Zgjt) 



m = -—^^l-z^,{t)sm5Mt) 



6(f)b{t) = AEb,a + ^ NbZb{t) + — NaZa{t) 
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2Kh 



Zb{t) 



^ V^^W) 



COS (506 (t) . 
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(22) 



where, 



AE, 



^,3 



E^L - EjR , UjiL - UjiR j^j- UjjL - UjjR jy 



E\ 



]a 



h 

dr 



2h 



2h 



h^ 



2m j 



\V^U^)\' + ^%V{v 



K, 



- / dr 



h^ 



2m. 



V<l>,x(r) • V<l>,^(r) + $,i(r) V{v) <l>,^(r) 



Uijo. = gij J dr ^l{r)^%{r) 



(23) 



with i,j = a,b and a = L,R. Let us consider a mixture with the same atomic mass for 
both components M = nia = m^t, which are trapped in the same symmetric double- well 
potential. Then, the localized modes are the same for both components but depend 
on the site: $i(i?) = $a,L(/?) = ^bMR)- Therefore, E^^ = E^,^ = E^j, = E^j, = E, 

UaaL = UbbL = UaaR = UbbR = U and UabL = UbaL = UabR = UbaR = U, Ka = Kb = K. 

Defining for each component the population imbalance and phase difference between 
both sides of the barrier, 



the above equations can be rewritten as: 

4(^) = -^R^\- Zl{t)'^\wb<\)a{t) 

Za{t) 



(24) 

(25) 



• , , NaUZait)+NbUZb{t) , 



Vi^4^) 



Zb{t) = -ujn\Jl- zl{t) sin (506 (t) 

■ N^Uz^it) + NailZait) 
d(pb{t) = h Ur- 



Zb{t) 



COS (50a (t) 



cos(50b(t) 



where ur = 2K/h is the Rabi frequency, the same for both species. It is useful to define, 
A = NU/huR, A = NJJ/huR, fa = Na/N, fi, = Nh/N and rescale the time as t -)■ ujRt, 

4(t) = -,/l-zlit)sm6Mt) (26) 



6Mt) = faAzait) + fbAZbit) + 



Za{t) 



v^T^^M 



Zbit) 



zl{t) sin (50b (t) 



5'<Pb{t) = fbAZhit) + fa~AZa{t) + 



Zb{t) 



COS S(f)a{t) 



cos(50b(t) . 



'1 - 4it) 

These equations correspond to two coupled nonrigid pendulums. The stability of these 
systems of equations have been analyzed recently in Ref . [23] . 
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3. 5. Improved two-mode model for the binary mixture 

As was noted for the scalar case in Ref . [9] , it is not mandatory to neglect any of the 
overlaps to obtain a closed set of equations relating the population imbalances and 
phase differences for a symmetric double-well potential. The complete set of two-mode 
equations were called the improved two-mode (I2M ) equations. 

In principle, if the experimental setup is appropriately chosen the left and right 
modes may be quite well localized at each side of the trap. In this case, the 
S2M equations are expected to provide quantitative agreement with the experimental 
data. When the two-modes are not so well localized, then it becomes necessary to 
consider the I2M to have quantitative agreement. In [9] the authors considered explicitly 
the set up of the Heidelberg group and showed that the I2M is necessary in the single 
component case to provide a quantitative understanding of the experimental data. 

Following similar steps as in the previous section and assuming the double-well 
potential to be symmetric as in the experiment, then the wave functions for the ground 
state and first excited state, $j±(r), have a well defined parity. The symmetry properties 
and the ortho-normalization conditions are capital to derive the coupled equations within 
the I2M model: $j±(r) = ±$j±(-r) , {^ia\^jp) = 5ij Sa/s , fori,j = a,b and a, /3 = 
+, — . The I2M provides an exact description of the dynamics in the symmetric double- 
well potential, with no approximations beyond the assumption of a two-mode ansatz of 
the total wave function \E'j(r;t), Eq. fl20|) . 

The resulting system of equations relating the population imbalance and phase 
difference for each component within the I2M approximation readdj]]: 



2K, 



ah 



Za{t) = ^y^ ~ ^a(^) sin(50a(t) 

5<Pa [t) = —T- + ^ /, .,,. COS 5<i)a{t) 



hit) = -^^^^l-zi{t)sin6Mt) 



h 

5m =^ + -^-^====cos6Mt) (27) 



with 



A,(t) = 2 77_iV, za{t) + 2 77J'^+„ N,Zh{t) , 

Afe(t) = 2 7j_ N,z,{t) + 2 7^^!"+_ N^zait) , (28) 

where we have defined 

1% =9^, /rfr$L(r)<l>5,(r) (29) 



aabb bbaa 

'+-+- ~ ' + -+ 



gab / rfr $„+ (r ) $„_ (r) $;,+ (r ) $5_ (r) , 



II Our system of equations differs slightly with the previously derived ones, of. appendix of Ref. |25) . 
We believe their system has some minor errors, which do not affect their discussion which is based on 
the S2M equations. 
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and 

2Kab{t) = (/i^ -l^\) + \ Na (7++ - 7--) (30) 

+ N, (7^'+ - 7-- - it- + 7-+) 

- Na[il\ + 7-- - 2^l^_)^l-zl{t)cos5Ut) 

- N,[^t+ + 7-- - it- - l%)\Jl-zl{t)cos5Ut) 

ix'j^ and ^■'_ are the chemical potentials of the ground and first excited state of the j 
component, that can be calculated from the time-independent GP equation for $j±, 
respectively. Analogously one can define 2Kba by exchanging the subindex a and h in 
the previous expression. 

Notice that we have kept the full 3D dependence of the wave functions <l>j±(r), 
instead of averaging the transverse spatial dependence as in Refs. [9l [25]. Thus, the 
coupling parameters gij in Eqs. ( 129|) are the 3D ones and are not renormalized. 

The equations for the I2M are essentially similar to the S2M . The main difference 
is that the tunnehng term, Kab(t), is time dependent and contains effects due to the 
interactions. As expected, if the localization of the modes is increased, i.e. by increasing 
the barrier height, Kabif) approaches the constant value, 2Kah — > /i" — l^\i which equals 
2K of Eq. ( I23l) . The coupled equations obtained in the I2M model reduce to well-known 
dynamical equations in two limiting cases: 

i) Setting to zero the overlapping integrals that involve mixed products of left and 
right modes of order larger than 1, the I2M equations reduce to the S2M model 
for the mixture, Eqs. (!25|l . 

ii) Assuming a noninteracting mixture, the inter-species interaction is Qab = 0, and the 
I2M equations for the mixture reduce to a two non-coupled system of equations, 
that are the dynamical equations of the I2Mfor a single component. Sec. 13.21 

As discussed in the introduction we are interested in the particular case of a 
binary mixture made of atoms populating two different hyperfine states. Then, both 
components have the same mass M, and are trapped in the same symmetric double-well 
potential. We initially restrict to the case in which the inter-species interaction is also 
almost equal to the intra-species one, g = Qaa = Qbh ~ Qab- This is the situation for 
F = 1 m = ±1 of ^'^Rb. This case allows straightforward comparisons between the 
results of the I2M and the ones obtained by solving the NPSEor GPID for a mixture 
explained in Sec. HI 

The ground and first excited states in a symmetric double-well potential are the 
same for both components. Moreover, since g = Qab the overlap integrals f l29|) reduce to: 



ll\ 


=1% 


= it^ - 7+^ 


iT. 


=lt'^ 


_ afc _ 


iT- 


=lt- 


^ ^ab ^ ab 



,aabb 



1+- , (31) 
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and the chemical potentials jj^^ = jj'^ = n^ with a = +, — . This yields the following 
relations: Kab = K^a and A^ = A;,. The I2M system reduces to: 

■ 2{NaZa{t) + Ni,Zk{t)h+. 2Kab{t) Zait) 

^ ^ ^/^-Zl{t) 



Zbit) = -^^|^^l-4(t)sin<50, 

.• , ,.. 2{N,Za{t)+NbZb{t))^+. 2Kab{t) Zbjt) 

0<Pb\t) = r + r , „, , COS d(pb{t) . 32 

h h ^i-z',{t) 

In this case both components obey the same system of coupled differential equations. 
Then, if the initial conditions are the same for both, Za{0) = Zb{0) and (50a(O) = (50(,(O), 
they will evolve with the same imbalance and phase, and no mixture effects will be 
observed. 

3. 6. Regimes for binary mixtures 

We proceed now to analyze the stability of the system of equations fl26l) . cf. see the 
appendix of Ref. [22] • As in the single component case, and in order to get analytical 
results that allow for a physical insight, we perform the study in the framework of the 
S2M approximation. First we note that an stationary point, defined by the equations: 
ij = and (50^ = 0, necessarily fulfills, 

sin 5(t)a = =^ 50° = 0, TT 

sin (506 = ^50|^ = O,7r, (33) 

and the following system of equations, 

1 
(1^ cos 500 



z^i 





zl = -zl i + -. ^ ^ 1 . (34) 

Therefore there are four different cases: (50° = 0,50° = 0), (50° = 0,50° = tt), 
(50° = TT, 50° = 0), (50° = 71, dcp^ = vr), noting that in all of them there is an obvious 
stationary point, z^ = z^ = 0. These stationary points will be referred to as "trivial 
stationary points" . We need to find the conditions for non-trivial solutions in each case. 
The stability of the system is analyzed by considering small variations around the 
stationary points for each of the four situations. Defining the displacements rji, 

Zait) = zl + riait) 

Zb{t) = zl + r]b{t) , (35) 

the following system of equations for the r/'s can be derived from Eqs. (1261) 



"^^ ) = -fi2 ( ^» ) (36) 
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l^cPl) ul/ul 



ul/ul 



(0,0) 

(7r,7r) 

(0,7r) 

(vr.O) 



where, 



1 + f l + V(/--/^)' + 4/a/6(A/A)2 



1 + f 1 



(/a-/.)2 + 4/,/,(A/A)^ 



1- 



l-V(/a-A)2 + 4/,/,(A/A)2 



1- 



l + V(/a-/.)2 + 4/,/,(A/A)^ 



1 + f U-/a) + Vl-4/./,(A/A)2 



1 + f U-/a)-Vl-4/./,(A/A)2 



1 + f (/a-/6) + \/l-4/./,(A/A) 



1 + f (/a-/6)-Vl-4/,/,(A/A) 



Table 2. Square of the frequencies of the eigenmodes of the S2M system, Eqs. (|26l) . 
Hnearized around the trivial stationary points, z'f — Q for the four different 5(j)^ 
combinations. 



((50a, 5(t)b 



(0,0) 

(7r,7r) 
(0,vr) 
(7r,0) 



ul/ul 



1 + A(l + 2/3/,/,) 
1 + 2A/3/J, 



1 



fa-fb 



l + {fa- ft)A 



fa—fb 



ul/ul 



1 - 2KPUh 

1 - A(l + 2/3/„/,) 

l + (/,-/„)A + ^g^ 



1 + 



fa-fb 



Table 3. Same as Table[2]but retaining up to the first order in /3, where A = A(l + /3). 
We assume fa > ft- 



n' 



OJf 



+ ^l 



1 + ifaAz^, + fM? 

1 + Ua Azl + f,K zlf 

fA Vl-{z'ay C0s6<Pl faAy /l-{z^,y COs6 

f,Ay/l-{z^,ycos6<p^, ft,A^l-{z^,ycos6. 



(37) 



In Tabled we give the exphcit values of the eigenfrequencies of Q for the trivial stationary 
points, zf = 0. These are obtained for the case under consideration, where A > and 
A>0. 

Approximate simpler expressions for the same eigenfrequencies can be derived for 
the case when A ~ A. Defining A = A(l + /3) and retaining up to terms of order /3, one 
obtains the frequencies listed in Table |3l 



3.6.1. Stationary points with (50° = 0, 50° = 0) 

In this case, the condition for the existence of non-trivial solutions to the equations 
depends on the slope at the origin of the two curves ([31]) [22] • The condition 

'A 1 \ ^A 1 , 

- + ^ 1 <1 



,A f,Aj VA faAj 
guarantees the existence of two additional solutions besides the trivial one. If we restrict 
ourselves to the case of A ~ A > 0, we have that ( l38l) cannot be fulfilled and therefore 
the only stationary point is the trivial one, z° = z° = 0. 
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In this case, it is straightforward to substitute the stationary point, z^ = z^ = Q 
and 5<pl = 5<j)l = Q into Eq. ([37]) to get, 

which has two eigenvalues, hsted in Table El 

In the very polarized case, fa ~ 1,/?, ~ 0, the population imbalance of the most 
populated component decouples from the less populated one and oscillates with the 
Josephson frequency wj = ui. The less populated component is driven by the other 
component and follows its dynamics, thus giving rise to "anti- Josephson" oscillations. 
The smaller frequency oscillation seen in the population imbalance of the less populated 
component is U2, which in this case is very similar to uji [30] . 

Also interesting is the non-polarized case, fa = fb = 1/2) then (assuming A ~ A, 
which is the case for ^^Rb), 

4(t) = - A/2{Za{t) + Zb{t)) - Za{t) , (40) 

Zi,{t) = - A/2{Za{t) + Zh{t)) - Zb{t) . 

and defining Az{t) = Zait) + z^it), 5z{t) = Za(t) — Zb{t) we have, 

Az{t) = -(A + l)Az{t) , 6z{t) = -6z{t) . 

Therefore, Az behaves as the single component case, oscillating with the usual Josephson 
frequency, wj = corV^ + A while 6z oscillates with the Rabi frequency, as would a single 
component case in the absence of atom-atom interactions. This mode can be further 
enhanced by imposing that Za{0) = —Zh{0) thus forcing both imbalances to oscillate 
with the same frequency. 

We have proposed in Ref. [30] to use these two configurations to extract the 
frequencies governing the dynamics of the system in order to obtain the microscopic 
atom-atom interaction. The idea was to profit from the fact that the difference between 
the inter- and intra-species interaction is small for the case of '^^Rb, A = A(l + /3), 
so we can use the expressions listed in Table El w^ = w|.(l + A(l + 2f3fafb)), and 
U2 = w^(l — 2Af3fafb)- Note that in the anti- Josephson case the oscillation with 
larger period is uj = w|(l + 0{Pfb)) and the shorter is uf = a;|j(l + A + 0{Pfb)), 
with (3 « 1 and fh « 1, allowing to extract both the Rabi and Josephson 
frequencies with good precision. The second configuration only has one frequency which 
is uf = uj]i{l + A(l + (3/2)) which allows to isolate the value of (3. 

3.6.2. Stationary points with ((50° = vr, (50° = 0) 

In this case, the condition for the existence of three stationary points is, 

,A M/ VA f,A) 

For the case considered here, A ~ A, and, in most applications, A > 1. Therefore, an 
appropriate choice of fa can ensure the existence of three stable points. The stability of 
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the trivial solution is checked by studying, 

,2 2 /I 0\ 2 ^ fa^ faA 



whose eigenvalues are listed in Table [2l The stability of the other two solutions is easy 
to study with the same tools. Simple analytic expressions are only attainable for the 
case A = A. Then we have, 

,_ , fl+K{faZl + hzlf 

1 i — UJ I 



R \ n 1 I A ('/ -vO I -f ^o\2 



1 + KUaZl + hzlf 



whose eigenvalues are, 

C^l = U^U^VaZl + hzlf) , i^l = cola + ^\faZl + hzlf) • (44) 

3.6.3. Stationary points with f(50° = vr, 50° = vrj 

The condition for the existence of three stationary points is in this case [22] , 

'A 1 \ /A 1 \ 

< 1 . (45) 



A m; va m. 

The eigenvalues corresponding to small oscillations around the trivial point are listed in 
Table El Its dynamical stability depends on the specific values of /j. A, A and ujr. For 
the case A = A, it is stable provided that u^ > A. 

The eigenfrequencies for the non-trivial solution are the same as for the case 
{6(f)a = 0, (506 = vr). For the simplest case, A = A, they are, 

Cof = CoUAVaZ'a + hzlf) , i^l = U^U^ + ^^faZ^ + hztf) . (46) 

4. Effective ID mean field approaches 

In the experimental realization [1] the condensate is confined by an asymmetric harmonic 
trap, characterized by Ux^ujy, and Uz, with a barrier on the x direction. Thus, in a 
first approximation one can assume that the dynamics takes place mostly along the x 
axis and derive descriptions of the system where the other two dimensions have been 
integrated out reducing the GP3D equation to an effective ID equation. There are 
different procedures to derive effective one dimensional GP-like equations starting from 
the three dimensional one. Their generalization to binary mixtures, with two coupled 
GP equations, or spinor BEG, with three or more coupled GP equations, is presented 
below together with the single component case. 

4.1. One dimensional Gross-Pitaevskii-like equations (GVIY) ) 

Assuming that most of the dynamics occurs in the direction which contains the barrier, 
the X direction in our case, one can approximate the wave function of the system by 

^(x, y, z- 1) ~ *i°(x; t) ipg.sXy) ^g.sXz) , (47) 
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where ^Pg,s. are the corresponding ground state wave functions for the trapping potential 
in the y or z direction in absence of interactions (in the case of harmonic traps they are 
Gaussian). In this way it can be shown [35] that \l/^^(a;; t) fulfills a Gross-Pitaevskii-like 
ID equation, 

d^^^{x;t) 



ih- 



dt 



2m 



^^^(x;t) 



(48) 



where the corresponding ID coupling constant is obtained rescaling the 3D one, giD = 
g/{27ra']_), with a± the transverse oscillator length, a± = ^/filmuj^, with uj^ = ^UzUJy- 
The extension to binary mixtures (and also to spinor condensates [36]) may be 
written down readily. 



ih 



dt 



-—dl 
2m ^ 



V{x)+J29ar,iBN,\^]''{x;t)\' 



j=a,b 



^l^{x;t) 



'^ dt 



-^''^^^^^ 



Y,9br,iBN,\^>f{x-M 

j=a,b 



^l^{x;t) 



(49) 



where, the rescaled couplings are gij-iB = gij/{2Tia] 



4-2. Non-polynomial Schrodinger equation (INPSE^ 

A more sophisticated reduction that includes to some extent the transverse motion 
of the elongated BEG in the corresponding potential is the so-called non-polynomial 
Schrodinger equation, proposed for a scalar BEG in Ref. [21] ■ The NPSE recovers the 
previously discussed ID reduction in the weakly interacting limit, but it has been shown 
to provide the best agreement with the experimental results on Josephson oscillations 
between two coupled BEGs [20]. The NPSE for the scalar case reads, 

.a^(x;t) 



ih- 



dt 



2m " ^ ' v^l + 2a,A^|^(x;t)|^ 



+ 



hjJx_ 



+ ^l + 2asN\^{x-t)\^ 



(50) 



^(x;t) . 



2 \^^l + 2a,A^|^(x;t)| 
The generalization of the NPSE for two components in a binary mixture of BEGs has 
been addressed in Ref. [37]. The system of equations, which become rather involved, can 
be greatly simplified in the case when all the interactions, both intra- and inter-species, 
are equal: 



ih 



d^!j{x]t) 
'di 



h^ 



2mi 



dl + V + g^D 



p{x]t) 



a/1 + 2asp{x;t) 



+ 



hu_i 



1 + 2asp(x; t) 



^,(x;t) (51) 



2 \ ^l + 2asp{x-t) 
where p{x;t) = A^a|\l'a(x; t)p -|- A'fe|\E'f,(x; t)p, j = a,b, and, as before, gm = g/{2Txa\), 



g = gaa = Qbb = gab = gba, and / (ix|^j 



.X 



1. 
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Figure 5. (left) Depiction of the potential in the x direction in units of h. The 
horizontal lines correspond to the single particle eigenenergies of the single particle 
Haniiltonian. (right) The first four single particle modes corresponding to the energies 
depicted in the left. 



5. Numerical solutions of the 3D Gross-Pitaevskii equation: single 
component 

Before analyzing the binary mixtures in the next section, we will present here numerical 
results for the single component to illustrate the main differences between the various 
two-mode models and ID reductions. 

As discussed in the introduction, we consider the same setup and the same trap 
parameters as in the experiments of the Heidelberg group |1]. There, a condensate 
of ^^Rb with 1150 atoms is confined to a fairly small region of ~ 5/im through the 
potential, 

V{v) = \m{uIx' + uy + uy) + V, cos\nx/qo) (52) 

with Ux = 27T X 78 Hz, Uy = 27t x 66 Hz, Uz = 27f x 90 Hz, go = 5.2/xm, and Vq = 413 h 
Hz. In Fig. |5]we show the potential in the x direction together with the first four energy 
levels of the single particle Hamiltonian and the corresponding modes. The energy 
levels of the single particle Hamiltonian show a clear separation between the two first 
eigenvalues, ground and first excited state, which are almost degenerate, and the next 
two. 

The atom-atom interaction strength is in this case, g = Anffa/M. The scattering 
length for ^^Rb is a = lOO.STas, therefore g/h=0M878 KHz fim^. Noting that the 
number of atoms is known up to 10% in the experiment, the relevant product, gN/h is 
in the range [51.22, 60.98] KHz fim^. Ref. [9] uses a value of 58.8 KHz fim^ to simulate 
the experimental setup. This large value of gN corresponds to a situation similar to 
panel (c) of Fig. HI where the possible dynamical situations we can have are: Josephson 
oscillations, i.e. closed orbits around the stationary point (2;°, 50°) = (0,0), and self 
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trapping regimes, usually funning phase modes. 

In the experiments, the system is prepared in a slightly uneven double- well potential 
which produces an initial population imbalance between both sides of the barrier. At 
t = the asymmetry is removed and the BEC is left to evolve in a symmetric double- 
well potential. In our numerical simulations the initial states with either (5(/)(0) = or 
TT are constructed in a different way than in the experiment. We build initial states 
which are by construction two- mode-like. First, we obtain numerically the ground and 
first excited states of the condensate in the double-well potential by solving the time 
independent GP equation (both for the ID reductions and the 3D case), then use those 
to build the left and right modes, Eq. (jlj), and finally construct initial states of any 
given initial imbalance, Zq: \E'2(,(r;t = 0) = acpiir) + e'*'^/30ij(r), with a^ + P"^ = 1 and 
a^ — P"^ = zq. The ground and first excited states are obtained by a standard imaginary 
time evolution of the equation from an initial state with the proper parity. The density 
profiles of the ground, first excited and left and right modes computed numerically are 
plotted in Fig. [H As can be seen, the left /right modes are indeed well localized at each 
side of the barrier. 

From these ground and first excited states we compute all the parameters entering 
in the S2M and I2M descriptions presented in Sees. 13. ![ and 13.21 The actual values of 
the parameters are, K/h = 0.00799 KHz, and NU/H = 1.19841 KHz for the S2Mand 
A/h = 1.19372 KHz, B/h = 0.03683 KHz, and C/h = 0.0023590 KHz for the I2m3. 
The values of the overlaps are: N-f++/h = 0.581746 KHz, N-f+_/h = 0.59803 KHz 

and A^7 /h = 0.623769 KHz. These numbers are used to generate the comparisons to 

S2Mor I2Min the following figures. 

In the full GP3D simulations we define the number of atoms in the left well 
as: Nilt) = j_ dx j_ dyj_ dz |\E'(r;t)p. The number of atoms in the right well 
is computed as Nji(t) = N — NiiJ:). From these values, the population imbalance 
reads, z(t) = (Niit) — Nfi{t))/N. Analogous definitions are used in the GPlDand 
NPSE equations. 

The phase difference between both sides of the potential barrier is computed in the 
following way. The phase at each point at a certain time, (f){x,y, z]t), is: 

^(x, y, z; t) = ^/p{x~y^z]¥) eyi.Y>{i 0(a;, y, z\ t)) , (53) 

where the local density, p(x, y, z; t) = \'^{x,y,z;t)\'^ . 

Averaged densities are defined as, i.e. integrating over the z component, 

/oo 
dz p{x,y,z]t) . (54) 

-oo 

To visualize the phase coherence along some of the planes we define, e.g. integrating 
the z component, 

1 r°^ 

(t){x,y]t) = — -/ dz p{x,y,z;t) ^{x,y,z;t). (55) 

p{x,y]t) J_^ 

% These values compare reasonably well with the ones provided in page 33 of Albiez PhD thesis [ID] , 
there they are given in units of ujx: Ajuj^ = 2.43572, Bjuj^^ ^ 0.0751497, CjuJx = 0.0048, and 
Kjujx = 0.0163. 
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Figure 6. The two smaller plots above depict in solid-black line the GP3D time 
evolution of z (left) and (50 (right), computed as explained in the text compared to 
the I2M predictions in dashed-red. Then we show 3D pictures complemented with 
contour plots, left, oi p{x,y]t), p{x,z;t) and p{y,z]t) at three different times, 0.5 ms 
(upper), 30 ms (middle) and 75 ms (lower), respectively. On the right of each plot we 
present a contour plot of the averaged quantum phase 0(x, ?/; i), 0(x, z; t) and <j}{y^ z; t) 
at the same times. They correspond to the first run presented in Fig. [SJ z(0) = 0.1 
and (5(/>(0) = 0. 



The phase on the left, 0L(t), is defined as, 

1 



/O foo 



dy / dz p{x,y,z;t) (j){x,y, z;t) 



(56) 



The phase on the right is defined accordingly. 

The way to implement the above averages over the phase has been done in the 
following way, 

iZo ^^ Im[^(x, y, z; t)] p{x, y, z; t) 



(j){x, y; t) = arctan 
(piif) = arctan 



/!! dz Re[^(a;, y, z\ t)\ p{x, y, z; t) ' 

/-oo dx /^ dy /^ dz Im[^(x, y, z; t)] p{x, y, z; t) 

lloo dx JZo dy IZo ^^ Re[^(x, y, z; t)] p{x, y, z; t) 



(57) 
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Figure 7. Similar to Fig.[6]but for a self-trapped case, z(0) — 0.7, (50(0) ~ 0, for three 
different times, 10, 30 and 75 ms and showing the averages over z. We plot p{x, y; t) and 
contour plots. On the right panels we present contour plots of the averaged quantum 
phase, (/>(a:, y; t). The phase coherence of the condensates at each side of the barrier is 
clearly seen. 



5.1. GP3B results 

In Figs. Eland [7] we present full GP3D simulations for a Josephson regime and a running 
phase mode self-trapped case, respectively. These two figures clearly show two relevant 
aspects of the problem. First, it is clear that during the full time evolution, which 
covers up to t = 80 ms in the figure, the system remains mostly localized on the two 
minima of the potential. Therefore, the density has a two-peaked structure over the 
considered time period. Secondly, the atoms in each of the two wells remain to a large 
extent in a coherent phase during all times. This can be seen from the uniform color, 
constant phase, at each side of the barrier in the right panels of the figures. These two 
characteristics of the time evolution of the 3D Gross-Pitaevskii equation support the 
use of two-mode approximations. 

The modulation of the density profiles on the transverse direction is seen to be 
small, with a mostly constant quantum phase in the region populated by the atoms. 
This indicates that the transverse dynamics can be integrated out to a large extent, as 
is done in the ID reductions discussed in Sec. HI 

The Josephson dynamics. Fig. El is clearly seen in the small upper panels depicting 
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Figure 8. Dynamical evolution of the population imbalance, z, between both sides 
of the barrier for a single component condensate. Solid (red) line corresponds to the 
GP3D , the dashed (blue) line to the NPSE , and the dotted (black) stands for the 
GPID. Panel (a) contains 50(0) = cases, with z(0) = 0.1,0.35, and 0.6. (b) 
Corresponds to the critical value, z(0) = 0.39 and 54'{0) = 0. (c) Depicts two self- 
trapped cases with an initial 5<j){Q) = tt, with z{Q) = 0.2, and 0.4. 



z{t) and 6(j){t). They both oscillate with the same period but with a phase-shift of 7r/2. 
A self-trapped case is shown in Fig. [71 The atoms remain trapped mostly on the 
left side of the trap (they start with an imbalance of z{0) = 0.7) and remain trapped 
in this potential-well during the considered time evolution. The coherence of the phase 
at each side of the potential barrier can also be appreciated in the figure, although here 
we should note that the right side of the barrier, being less populated, is concentrated 
on a smaller {x, y) domain. 

5.2. Comparison between the different models 

The GP3D cases described above indicate that within the configuration considered 
here the two commonly employed two-mode models and ID equations are expected to 
be reasonable. In this section we present comparisons between the different approaches 
described in the previous sections: ID reductions (NPSE , GPID ) and two-mode 
models, S2M and I2M . 



5.2.1. GP3D vs ID reductions: GPID and NPSE 

In Fig. IHlwe present the time evolution of the population imbalance for the different 

Josephson, and self-trapping. 



dynamical conditions described in Sec. 13. 3[ i.e. Josephson, and self-trapping. We 
compare the full GP3D (solid red) with the two previously described ID reductions, 
GPID (dotted black) and NPSE (dashed blue). 

First, we note that the dynamics emerging from the GP3D is indeed similar to 
what was predicted by analyzing the S2M equations in Sec. 13.31 Qualitatively, the 
GP3D simulations do follow the patterns predicted by the two-mode approximations. 
Lets us briefly describe each of the results: 

a) The first panel, (a), contains simulations performed with zero initial phase 
difference, i.e. Josephson oscillations and self-trapping cases. For the Josephson 



Weakly linked binary mixtures of F = 1 Rb Bose-Einstein condensates 



28 



t=50 ms 



0.4 

0.2 


0.4 

0.2 



z(0)=0.1,5(^(0)=0 




-L 



z(0)=0.386, 5(^(0)=0 




z(0)=0.35, 5(^(0)=0 



!_. 




z(0)=0.2, 5(|)(0)=7t 




z(0)=0.6, 5(^(0)=0 



1^ 




_L 



z(0)=0.4, 5(^(0)=7t 




-3 3 
x()xm) 



-3 3 
X ()xm) 



6 -6 



-3 3 
x()am) 



Figure 9. Snapshots of the axial density profiles, p{x;t) (fim)^^ at t = 50 ms 
calculated by means of the GP3D evolution (solid red line), the NPSE (dashed blue 
line), and the GPID (dotted black line). The initial conditions correspond to the ones 
used to generate Fig. [Sj 



cases, z{0) = 0.1,0.35, the imbalance oscillates with a frequency which is mostly 
independent of the initial imbalance (for small imbalances). With z{0) = 0.1 the 
oscillations are almost sinusoidal, while as we increase the initial imbalance their 
shape becomes more involved but remaining periodic. In the self-trapped case, 
z{0) = 0.6, the atoms remain mostly on the initial side of the trap and there are 
short and small periodic oscillations as predicted by the two-mode models. At 
longer times, the imbalance is seen to decrease smoothly, implying a departure 
from the predicted two-mode dynamics |38] . 

The two ID reductions give qualitatively similar results in most situations to 
GP3D , but not quantitatively in all cases. The NPSE is seen to reproduce 
very well the GP3D in all the runs up to times near ~ 40 ms. Above those 
times, the period of oscillation predicted by the NPSE is slightly shorter than the 
GP3D one. The GPID on the contrary only captures the amplitude of oscillation 
in the Josephson cases, failing in all cases to give the same period as the GP3D or 
the NPSE. Moreover, the GPID departs notably from two-mode for the self- 
trapped case. It does predict self trapping, but more than two modes contribute to 
the time evolution. 

b) Panel (b) is computed near the critical value of the full GP3D , z{0) = 0.39 for 
50(0) = 0. The GPID and NPSE predict a critical initial imbalance close to the 
value predicted by the GP3D . 

c) Panel (c) contains two self trapped cases obtained with an initial 50(0) = vr and 
^(0) = 0.2, and 0.4. Notice that for 50(0) = it the critical imbalance is smaller. 
The discussion is similar to the Josephson case, i.e. the NPSE captures most of 
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Figure 10. Dynamical evolution of the population imbalance between the two sides of 
the barrier for a single component condensate. The GP3D (solid red) is compared to 
the I2M (dashed blue) and the S2M (dotted black) results. The parameters entering in 
the two-mode descriptions are given in the text. Panel (a) contains runs for (50(0) = 0, 
with z(0) = 0.1,0.35, and 0.6. (b) Corresponds to the critical value for z(0) = 0.39 
and (50(0) = 0. (c) Depicts two self-trapped states obtained by an initial (50(0) = tt, 
with z(0) = 0.2, and 0.4. 



the dynamical features of the GP3D while the GPlDonly provides a qualitative 
understanding of the problem. 

These results justify the use of the NPSEin Ref. [4J to analyze their experiment. 

To further explore the quality of the ID reductions, we present in Fig. [9] the 
density profiles in the x direction after integrating the y and z components, p(x; t) = 
/_ dy j_ dz |\l/(x, y, z] t) p at t = 50 ms. The agreement between the NPSE and the 
GP3D is very good in most situations, except for the critical case, as expected. In 
all cases the density profiles show a clear bi-modal structure. The GPID , as could 
be inferred from the previous results, does not predict the correct density profiles and, 
as seen in the self-trapped case, (-^(O) = 0.6,(50(0) = 0), do show the contribution of 
higher modes. The critical initial imbalance starting with no phase difference that we 
find numerically by means of the GP3D is the same as found in Ref. [9], Zc = 0.39, and 
differs from the one reported in Ref. |1], ^c = 0.5. 



5.2.2. GP3D ws two-mode approximations, S2M and I2M 

As explained above, the use of two-mode models is suggested by the GP3D results, 
see Figs. [6] and [71 What is, a priori, not clear, is whether the extra assumption used in 
deriving the S2M (which are the most commonly employed equations) will work for each 
specific double-well potential. As discussed in Sec. 13.21 the conditions of the Heidelberg 
experiment are such that the S2M predictions are not good. This does not mean that 
the dynamics is not two-mode but that the overlaps involving high powers of the two 
localized modes are not negligible as assumed in deriving the S2M equations. 

In Fig. do] we compare GP3D (sohd red), the S2M (dotted black) and the 
I2M (dashed blue) results using the parameters calculated microscopically from the 
ground and first excited state of the GP3D . Both two-mode schemes predict the same 
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phenomenology and thus quahtatively capture the dynamics of the system. At the 
quantitative level, however, the I2Mis clearly better. In the run with 2;(0) =0.1 and 
50(0) = 0. (panel (a)), both the S2Mand I2M predict a similar behavior with the 
correct amplitude and oscillation period close to the GP3D one. As the imbalance 
is increased, e.g. (P] considers z(0) = 0.28), the S2M fails to describe the correct 
period and predicts smaller amplitudes. This is analyzed in full detail in Ref. [9]. The 
critical initial imbalances determined by both two-mode approaches are smaller than 
the GPSDone. For the latter they predict a self-trapped case, see panel (b). Finally, 
for the self-trapped cases with (50(0) = tx (panel (c)) the I2Mgive similar oscillation 
amplitudes with shorter periods than the GP3D . The S2M fails both in reproducing 
the amplitudes and the periods. 

6. Numerical solutions of the 3D Gross-Pitaevskii equations: binary 
mixture 

As discussed in Sec. [2], one feasible way of experimentally prepare binary mixtures of 
BECs is to consider a number of atoms populating the m = ±1 Zeeman components 
of an ^^Rb F = 1 spinor. The experimental observation of Josephson tunneling 
phenomena by the Heidelberg group seems to be possibly extended to trap both Zeeman 
components [M] . In this case the two components of the mixture have the same mass, 
M = nia = nib, and equal intra-species interactions, gaa = Qbb = 9- With respect to the 
inter-species interaction we will consider the case of ®^Rb which implies Qab ~ 9- 

The mean field GP3D system of equations governing the dynamics of the three 
components of an F = 1 spinor BEC can be written as |39j . 



i^— ^ = [Us + C2(n±i + no - n^i)]ip±i + Ca^o^^i ' 

i^— — = [Us + C2(ni + n_i)]i)Q + C22^i^^tjj_i , (58) 

with Us = —h'^/{2M) V^ + V^ + con being the spin-independent part of the Hamiltonian. 
The density of the m-th component is given by nm(r) = |'0m(r)p, while n(r) = 
Sml^m('^)P is t^^ total density normalized to the total number of atoms A^. The 
couplings are cq = 47r^^(ao + 202) /(3M) and C2 = 47r/i^(a2 — ao)/(3M), where oq and a2 
are the scattering lengths describing binary elastic collisions in the channels of total spin 
and 2, respectively. Their values for ^^Rb are oq = 101. 80^ and 02 = 100. 4aB |40j . 
Since the spin-dependent coupling, C2, is much smaller than the spin-independent one, 
Co, and the total number of atoms that we will consider is relatively small A^ = 1150, the 
population transfer between the different components can be neglected |i30j. Therefore, 
in our calculation the number of atoms in each sublevel remains constant in time allowing 
to treat the system as a real binary mixture of components a and b. Comparing the 
system of Eqs. ([1]) and (l58ll the value of the couplings can be read off, 9aa = 9bb = Co + C2 
and 9ab = 9ba = Cq - C2. 

Once the total number of atoms is fixed we want to investigate the Josephson-like 
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Figure 11. (left) Values of the frequencies, w/wij, listed in Table |3] for the specific 
conditions considered in the numerical simulations as a function of the fraction of atoms 
in the a component, fa- The notation is as follows, uJi-a^, with i = 1, 2 and a, /3 = 0, n. 
(right) The conditions for the existence of the non-trivial equilibrium points given in 
Eqs. (|551 SU H5)) . upper panel, as function of /„ for the conditions described in the 
text. The lower panel contains the explicit equilibrium points zj^, z^ as a function of 
fa obtained by solving equations (1341) . Note that each equilibrium point has a trivial 



partner which is obtained by flipping the sign of z^, z^- 



dynamics for different number of atoms populating each component A^^ = faN and 
Nb = fbN and for different initial conditions Za{0), Zb{0), (50a(O) and 50fe(O). 

The values of A = NU/Hur and A = NU/hhjji are A = 74.278 and A = 74.968. 
With A/A = 0.99. These are obtained from the microscopic 3D parameters computed 
in the scalar case, with the same total number of particles, Sec. |5l This is reasonable for 
the case we are considering where gaa = Qbb ~ Qabi which implies that the ground state 
wave functions for the GP equations of the mixture do not depend on fa and fb for a 
fixed total number of particles. This would certainly not be the case if Qaa = 9bb 7^ 9ab^ 
in such case one would need to recompute the ground state wave functions for a and h 
for each value of /„. 

Following the discussion in Sec. 13.61 where the predictions of the S2Mwere 
discussed in detail, the system has the trivial equilibrium points, listed in Table [3] 
with /3 = 0.009. In Fig. [11] we show the values of the two eigenfrequencies for each 
of the trivial equihbrium points hsted in Table |3] for the specific conditions described 
above. The figure shows a number of important features about the stability of the trivial 
equilibrium points. First, the (z°, (5(/)°, z°, 50°) = (0, 0, 0, 0) is always stable regardless of 
the total polarization of the system (measured by fb—fa)- Second, the (2;°, (5(/)°, z°, 50°) = 
(0, TT, 0, tt) mode is always unstable, as seen by the negative value taken by the square 
of the frequencies. Third, the (2;°, 50°, 2;°, 50^) = (0, 0, 0, vr) mode should be stable for 
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fa < 0.43, correspondingly the (2°, 50°, z^, 50°) = (0, vr, 0, 0) is stable for /^ < 0.43 and 
therefore there is a range of polarizations, given by 0.43 ^ fa ^ 0.57 where the only 
trivial mode which is stable is the (z°, 50°, z^, 54>i) = (0, 0, 0, 0). 

The non-trivial equilibrium points in this case can be obtained by analyzing the 
conditions given in Sec. 13.61 For (50°, 50°) = (0, 0) there are no equilibrium points 
apart from the trivial one, due to A ~ A. In the other three cases there are non- 
trivial equilibrium points depending on the specific values of fa- In Fig. [TlT right) we 
analyze their existence. First, we note that there are non-trivial points corresponding 
to (50°, 50°) = (0, n) provided fa < 0.37, correspondingly there are also equilibrium 
points for (50°, 50°) = (0, tt) if ff, < 0.37. There is also a non-trivial equilibrium point 
for (50°, 50°) = (tt, tt) regardless of /„. As can be seen in the figure, all these non- 
trivial equilibrium points correspond to fairly imbalanced conditions and can in most 
cases be understood in simple terms from the analysis of the scalar case. For instance, 
the equilibrium point for (50°, 50°) = (tt, tt) corresponds to 2;° ~ 2° ~ 1 (or —1), 
which can be understood as having both components locked in a vr-mode. Similarly, 
the equilibrium points in the (0, vr) or (tt, 0) cases exist whenever the most abundant 
component is populated enough to drive the dynamics close to being vr locked. 

6.1. GP3D calculations: phase coherence and localization 

The numerical solutions of the GP3D presented in Sec. |5]for the single component case 
showed two features. First, the atoms remained mostly localized in the two minima 
of the potential well and secondly, each group of atoms had to a large extent the same 
quantum phase. This, clearly supported the picture of having two BEC, one at each side 
of the barrier, with a well defined phase at each side during the dynamical evolution. 
Essentially those are the premises used to derive the two-mode models, both for single 
component and for binary mixtures, as we did in Sec. [3l 

As in the scalar case, our exact GP3D numerical solutions of the dynamics of 
the binary mixture in several initial conditions of population imbalances and phase 
differences show two distinctive features, see Fig. [121 First, the density of atoms for 
each component is always bi-modal, with the two atom bunches centered around the 
minima of the potential well. Secondly, the phase of the wave function is mostly constant 
for each species at each side of the potential trap. Thus, we find that the GP3D does 
predict the dynamics to be mostly bi-modal also for the binary mixture case. 

At the end of the section we will consider some deviations from the bi-modal 
behavior that are found in very specific conditions, e.g. for very large population 
imbalances and also when analyzing a case with Qab 7^ Qaa = dbb- 

6.2. Small oscillations around z°^ and 50°^ = 

The two predictions of the S2M described in Sec. 13.61 are confirmed by the NPSE and 
GPID simulations as can be seen in Figs. [13] and [Ml In Fig. [131 (left panels) we consider 
a very polarized case, /„ = 0.8. As expected from the two-mode analysis the dynamics of 



Weakly linked binary mixtures of F = 1 Rb Bose-Einstein condensates 



33 




X (|lm) 



3-113 



X {[J.m) 



Figure 12. Full GP3D calculations of the dynamics of a binary mixture with 
za(0) = 0.5, Zb{Q) = 0.2, ,50,(0) = 0, 5<Pb{Q) = 0, fa = 0.25 and fb = 0.75. The 
upper four plots depict, from left to right, Za{t), S(f>a(t), Zh{t), and d(j)b(t) in solid black 
compared to the I2M prediction, dashed-red . Then each row contains from left to 
right: 3D depictions complemented by contour plots of Pa{x, y\t), a contour plot of the 
averaged phase (l)a{x,y;t), 3D depictions complemented by contour plots oi ph{x,y]t), 
and a contour plot of the averaged phase 0b (x, y; t). Each row corresponds to a different 
time, .5 ms (upper), 20 ms (middle) and 60 nis (lower), respectively. 



the most populated component should to a large extent decouple from the less populated 
one and perform fast Josephson oscillations with a frequency close to the corresponding 
one for the scalar case, uj = ujRy/1 + A. The GP3D simulation is seen to confirm the 
above and follow closely the predictions of the I2M . The less abundant component is 
strongly driven by the most populated one and shows an anti- Josephson behavior as 
described in Ref. [30] . 

Another prediction is related to the behavior of Za + z^ and Za — zi, in the non- 
polarized case, fa 



fh. As explained in Sec. I3.6[ in this case the difference, Za — Zb, 
should enhance the long mode which oscillates with the Rabi frequency of the system, 
while the sum Za + Zb should mostly oscillate with the Josephson frequency. In the 
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Figure 13. Behavior of the population imbalance, Za{t) (solid lines), and Zb{t) dashed 
lines, and phase difference, 54)a{t) (solid lines) and 6(j)b{t) (dashed lines), computed 
using GP3D (black lines), NPSE (blue lines), and I2M (red lines) in a polarized case, 
fa = 0.8, left, and a zero polarization case, fa — 0.5, right, respectively. The initial 
conditions are Za{0) = 0.1, Zf,(0) = —0.15 and 6(j)a{0) = 6(j)b{0) = for the left panels, 
and Za{0) = —zi,{0) = 0.15 and S(j)a{0) — S<j)b{0) = for the right panels. 
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Figure 14. Behavior of the population imbalance in NPSE (red) and GPID (black) 
simulations in the zero magnetization case, fa = ft- The initial conditions are 
Za{0) — 0.1, Zb{0) = 0.2 and S(j){0) = 0. The upper panels correspond to (a) Za{t) 
(solid line) and Zb{t) (dashed hue) obtained with the GPID equations, (b) Za{t) and 
Zb{t) obtained with the NPSE equations, (c) behavior of Za{t) — Zb{t) for GPID (solid) 
and NPSE (dashed), and (d) behavior of Za{t) + Zb{t). 
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Figure 15. (left panels) Evolution of the population imbalance of each component for 
a binary mixture with fa = 0.25. The upper panel shows Za{t)^ and the bottom panel 
Zb(t). The solid (black) line corresponds to the I2M model and the dashed (red) line to 
the NPSE . The initial conditions are Za{Q) = 0.5, Zf,(0) = 0.2, 5(t>a{0) = 6^0) = 0. 
(right panels) As in the left panel, but with fa = 0.6 and initial conditions Zo(0) — 0.45, 
Zb{0) = -0.35, 6(j)a{0) = SMO) = 0. 



right part of Fig. [13] we present the extreme case when Za{0) = —-26(0) computed 
with GP3D , NPSE and I2M . In this case, both population imbalances and phase 
differences oscillate mostly with the Rabi frequency of the system, keeping during the 
time evolution Za + Zb ^ 0. 

As seen in Fig. [H] both ID reductions produce qualitatively similar physics. The 
only important difference is that the frequency of the Josephson oscillations is higher in 
the GPID , as occurred already for the single component, see Sec. |5l 

Interestingly, they predict different Josephson oscillations while the Rabi 
frequencies are similar. In panel (c) of Fig [H] the long oscillation corresponding to 
the Rabi mode is seen to agree well with the corresponding long oscillation seen in the 
right panels of Fig. [131 The Josephson-like oscillations of binary mixtures of spinor 
F = 1 s^Rb BECs around the (2°, 50°, z^, 50^) = (0, 0, 0, 0) are therefore essentially 
controlled by two frequencies, ur and uj. 

As a general statement, in the conditions of the Heidelberg experiment, as occurred 
for the scalar case, the I2M produces more reliable results than the S2M model, which 
are not shown in the figures. Notice that the parameters that we use for the I2Mare 
extracted from the GP3D calculation as given in Sec. [5] Other representative cases 
with (50a(O), 50fe(O)) = (0,0) but with larger initial imbalances, ^j(O) ~ 0.5 are shown 
in Fig. [151 On the left side of the figure we show the population imbalance of each 
component for a simulation with /„ = 0.25. In this case the dynamics is controlled by 
uj. The panel on the right depicts a simulation with fa = 0.6 and close to opposite 
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Figure 16. Two simulations with the same initial conditions, Za(0) = 0.1, Zb(0) = 
—0.15, (50a (0) — and (506(0) = tt but with different compositions of the mixture. 
The case on the left has fa — 0.2 while the case on the right fa = 0.8. The blue lines 
are obtained by means of a full GP3D , the black lines are the NPSE results, and 
the red lines are the I2M results. Solid and dashed lines correspond to the a and b 
components, respectively. 



initial population imbalances. In this case, both frequencies uj and ujr show up in the 
evolution. The I2M provides a satisfactory description of the dynamics. 



6.3. Small oscillations around z\, 



and (50° 



TC 



As explained above, for these conditions there can exist up to three stationary points 
depending on the specific value of /„ considered. The trivial equihbrium point exists 
provided fa < 0.43, see Fig. [11] This prediction of the two-mode models is observed 
both in the GP3D and NPSE as it can be seen in Fig. [161 In the figure, we consider 
a simulation with Za{0) = 0.1, 2fe(0) = —0.15, and /„ = 0.2 < 0.43 (left panels). The 
population imbalance (upper panel) of both components oscillates in the usual Josephson 
regime. At the same time, the phase difference oscillates with its characteristic phase- 
shisft of 7r/2 with respect to the imbalance (lower panel). The phase of the a component 
oscillates around 5(^a = while S4>b does oscillate around 6(f)b = tt. 

A completely different picture emerges when the fraction of atoms in both 
components is exchanged, fa = 0.8 > 0.43 (right panels), with most of the atoms 
populating the a component. In this case, the oscillation amplitude is large, both 
components remain trapped on their original sides and the phase difference becomes 
unbounded. This should be considered as a genuine effect of the binary mixture as each 
component follows a running phase mode at each side of the potential barrier. 
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Figure 17. Two simulations with the same initial conditions, Za{0) = —0.78, 
Zb{0) = 0.99, S(j)a — and S(j)b = tt, but with different composition. The case on 
the left has fa = 0.1 while the case on the right fa = 0.9. The red lines are obtained 
by means of a full GP3D while the black lines are the NPSE results. Solid and dashed 
lines correspond to the a and b components, respectively. 



The comparison between the NPSE and GP3D is very satisfactory. The 
NPSE captures almost completely the dynamics up to times of 100 ms. In all cases, 
the NPSE reproduces correctly both the phase difference and population imbalance. 
The only sizeable discrepancies occur for times > 80 ms in the run without equilibrium 
point (right panel). 

The I2M gives a good qualitative picture of both cases but fails to provide 
predictions as accurate as the NPSE , as happened in the scalar case, see for instance 
Figs. M and [lOl In particular the predicted periods of oscillation are much longer than 
the actual ones. 

An example of simulations around non-trivial equilibrium points is presented in 
Fig. [T71 As explained previously, these involve very large and opposite initial population 
imbalances for both components. In Fig. [17] we consider a case with initial conditions 
very close to the predicted equilibrium point using the standard two-mode, and described 
in Fig{TT|, Za{0) = —0.78, and Zh{0) = 0.99, with fa = 0.1. Also in the same figure we 
consider a similar run but with fa = 0.9. In both cases the NPSE and GP3D predict 
a very similar dynamics. These simulations will be discussed again in Sec. 16.51 as they 
exhibit effects which clearly go beyond a two-mode approximation. 



Weakly linked binary mixtures of F = 1 Rb Bose-Einstein condensates 



38 



1 


II 


0.5 

N 







7\f^y'^ " : 


-0.5 


— 




— 


80 


— 




— 


-e- 

CO 40 


^^^^ 





^--r^'^T^^^ 1 







50 

t (ms) 



100 




Figure 18. Two simulations corresponding to (left) Za(0) = 0.4, Zfc(O) — —0.2, 
,50a(O) = ^, (50b(O) = TT and /„ = 0.9, and (right) z„(0) = 0.9, Zb(0) = 0.85, 5<j)a{Q) = tt, 
i5(^h(0) = TT and /a = 0.9 The blue lines are obtained by means of a full GP3D while 
the black lines are the NPSE results. Solid and dashed lines correspond to the a and 
h components, respectively. 



6.Jf^. Small oscillations around z^^ and 50° ^ 
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The trivial equilibrium point is not stable in the considered conditions as seen in Fig. [TTJ 
The non-trivial one, however, is only attainable if extremely imbalanced configurations 
for both components are considered. This case would correspond essentially to having 
both components in a vr mode state, which in our conditions only exists for z ~ 1 as can 
be seen in the blue spots in panel (c) of Fig. HI In Fig. [18] we present two simulations 
with different initial conditions. First, we consider a simulation with Za{Q) = 0.4 and 
Zb{0) = —0.2, with fa = 0.9. The behavior is understood in simple terms, the most 
populated component remains self-trapped while the other component is forced by 
the other one. The phase evolves unbounded. The figure again contains GP3D and 
NPSE simulations. 

The second simulation (right panels) is closer to a non-trivial equilibrium point, 
we consider Za{0) = 0.9 and ^^(O) = 0.85 with fa = 0.9. In this case, both components 
remain self trapped, the phase difference is unbounded, but we do not get the expected 
behavior of two n modes because the initial imbalances are not close enough to z" ~ 1. 



6.5. Effects beyond two-mode 

Most of the dynamics described in the previous sections can to a large extent be 
understood within the two-mode models developed in Sec. [3l There are, however, a 
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number of situations where the two-mode fails. Some are a direct consequence of having 
two components evolving in the same double-well potential, others are due to having 
initial configurations, mostly with large initial imbalances, producing situations where 
the atom-atom interaction energy per atom is comparable to the gap between the first 
excited state and the second/third excited states. 

We can distinguish two different cases: (a) involving excitations along the 
coordinate which contains the barrier, (b) involving excitations of the transversal 
coordinates. 

An example of (a) is seen in Fig. [T71 There, as clearly seen in the density profiles 
along the x direction, the two-mode approximation is clearly not valid. The simplest 
way of seeing this is by noting the zero in the density of one of the components at 
X ~ 2/i m. This effect beyond two-mode is well taken care of by the NPSE which 
reproduces the density profile quite well during most of the time evolution considered 
in the simulation. Thus, the excitations of higher modes along the direction which has 
not been integrated out in the ID reduction do not pose a great difficulty to the ID 
reductions. 

The second type, (b), of effects beyond two-mode involve excitations of the 
transverse components. These effects are present in any binary mixture calculation 
whenever the intra- and inter-species interactions are not equal. To enhance this effect, 
and also to explore the interesting symmetry breaking phenomena described in Ref. [25], 
we consider a case with Qaa = dbb, but with gat = gta = '^■'^Qaa- Therefore, now the 
inter-species interaction strength is larger than the intra-species one. The two-mode 
prediction for this case, S2M , which was analyzed in Ref. [25] shows a large symmetry 
breaking pattern during the time evolution of the system. In Fig. [19] we consider a 
full GP3D simulation of a representative example with Za{Q) = —0.2, -2^(0) = 0.1, 
50,(0) = 50,(0) = 0, and /, = 0.7. 

The qualitative prediction of the I2M also shows the symmetry breaking, and the 
two components do separate from each other and mostly concentrate on one of the wells 
as time evolves. But, as it can be seen in the 3D depictions of p{x, y; t) at three different 
times, the evolution of the system departs almost from the beginning from the two- 
mode. At t = 1 ms we have the density distributions of each component corresponding 
to a small initial imbalance. Then at t = 11 ms, we can already see that the most 
populated component is expelling the other one from the minima of the potential. This 
fact can be appreciated as a four peaked distribution, pb{x,y;t). After that, each of 
the components start to accumulate on their original sides following qualitatively the 
prediction of the I2Mand thus presenting the symmetry breaking pattern discussed 
in Ref. [25]. The two- mode approximation is in this case broken for a short period of 
time, when the first modes along the transverse directions are excited due to the large 
inter-species interaction. 
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Figure 19. Full GP3D calculations of the dynamics of a binary mixture with 
z„(0) = -0.2, zb(0) = 0.1, <5(/.,(0) = 0, ,506(0) = 0, /, = 0.7 and h = 0.3. As 
explained in the text in this case, Qaa — Qbb and gab — '^■igaa- The upper two plots 
depict Za{t) (left) and Zf,{t) (right). Then each row contains from left to right: 3D 
depictions complemented by contour plots of Pa{x, y]t), a contour plot of the averaged 
phase (j)a{x,y]t), 3D depictions complemented by contour plots of pb{x,y;t), and a 
contour plot of the averaged phase (j)b{x,y;t). Each row corresponds to a different 
time, 1 ms (upper), 11 ms (middle) and 51 ms (lower), respectively. In all cases, solid 
black lines are computed with GP3D and dashed red ones with I2M . 



7. Conclusions 

We have presented a thorough investigation of the mean-field dynamics of a binary 
mixture of Bose-Einstein condensates trapped in a double-well potential. The rich 
dynamical regimes which take place in binary mixtures, like double self-trapped modes, 
Josephson oscillations, or zero and vr bound phase modes, have been scrutinized by 
performing full GP3D simulations covering all the relevant initial conditions. The 3D 
numerical solutions of the Gross-Pitaevskii equations have been used as a benchmark to 
critically discuss the validity of the most common ID reductions of the GP equations, 
GPlDand NPSE, and the often employed simple two-mode reductions, S2Mand 
I2M. 

The full 3D solutions of the binary mixture have shown to have a large amount 
of phase coherence and localization at each side of the potential barrier for both 
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components, predicting a dynamics which is mostly bi-modal. This feature permits 
to speak of Bose-Einstein condensates at each side of the barrier, where the atoms 
mostly share a common phase, and to support the use of two-mode approximations, 
which analytical solutions allow to gain physical insight into the problem. 

To fix the conditions of the dynamics, we have focused in one particular setup that 
corresponds to a natural extension of the experiments reported in Ref. |3]: the case 
of a binary mixture made by populating two of the Zeeman states of an F = 1 ^^Rb 
condensate. As discussed in the present paper, this setup already allows to observe and 
characterize a large variety of phenomena which are genuine of the binary mixture, e.g. 
anti-Josephson oscillations in highly polarized cases, long Rabi-like oscillating modes, 
zero and vr locked modes, etc. 

For the sake of completeness and to better frame the physics of the binary mixture 
we have provided a detailed description of the single component dynamics, with explicit 
expressions for all the commonly employed approximations to the 3D mean field Gross- 
Pitaveskii equation. The natural extension of the latter to the binary mixture, i.e. 
S2Mand I2M equations and ID reductions, have been consistently derived providing a 
self-contained reference, easy to read, with all the relevant formulae used in the article. 

The standard two-mode model, with its microscopic parameters computed with the 
GP3D , has been used to reexamine the existence and stability of the different regimes 
that can occur in both single component and binary mixture condensates, describing 
the Josephson oscillations and the macroscopic quantum self-trapping, including running 
phase modes and zero- and vr— modes. 

The comparisons between the two-mode models and the numerical solutions of 
the GP3D show an excellent agreement for conditions close to the stable stationary 
regimes predicted by the two-mode models. As we depart from those stable points, 
the S2M fails to provide a quantitative agreement with the results obtained with the 
GP3D equations. The range of validity of the I2M is much larger, fully capturing the 
dynamics of single and binary mixtures for a larger set of initial conditions. 

The two most commonly employed dimensional reductions of the GP3D , the 
GPID and NPSE , have been shown to differ substantially among each other, with 
the NPSE being clearly in much better agreement with the original 3D dynamics in a 
broader set of conditions. In general, the GPID describes essentially the correct physics 
but quantitatively far from the GP3D predictions. Also, for self-trapped cases already 
in the single component case, it departs from the two-mode behavior earlier than the 
GP3D or the NPSE. The agreement between the NPSE and the full 3D dynamics is 
astonishingly good both for single component and the considered binary mixtures, where 
the intra- and inter-species are very similar and the NPSE equations are particularly 
easy to handle. This agreement is not only seen on fully integrated magnitudes, for 
instance population imbalances, but also on the density profiles predicted along the 
direction hosting the barrier. 

We have also considered two situations where the two-mode approximation fails. 
This is naturally due to the excitation of higher modes. Two different cases have been 
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described, first the excitation of modes in tlie direction of the barrier and secondly, 
excitation of modes in the transverse direction. The NPSEhas been shown to capture 
perfectly the excitations along the barrier direction, reproducing the integrated density 
profiles obtained with the GP3D . The second case has been studied in a simulation 
performed with different intra- and inter-species, which can be achieved in principle 
experimentally through Feshbach resonance modulation of the scattering lengths. In 
this case, the dynamics of the less populated component in each side of the trap departs 
notably from the two-mode with clear excitations of transverse modes, seen already in 
the density profiles along a transverse direction. 

The present article is intended both to motivate the experimental effort to study 
binary mixtures of BECs, where we have shown that a large variety of phenomena 
related to phase coherence and localization can be observed, and to serve as a tool in 
the analysis of such experiments providing a concise and self-contained derivation of the 
most commonly used models. 
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